% TODO: update this.
%run this script after the full analysis has been completed. Load the
%analysis summary text file and plot all hydrodynamic fit parameters.

%this will need to be modified...

fname = 'ones_diffusion_summary.txt';
save_results = 1;

%read summary file
%titles
fid = fopen(fname, 'r');
titles = split(fgetl(fid), '\t');
fclose(fid);

%data
dat_mat = dlmread(fname, '\t', 1, 0);

Ts = dat_mat(:,1);
TsUnc = dat_mat(:,2);
Dens =  dat_mat(:,4);
DensUnc = dat_mat(:,5);

%hydro parametrization #1 (Gamma, Dmicro, c, nu)
Gammas = dat_mat(:,7);
GammasUnc = dat_mat(:,8);
Deffs = dat_mat(:,9);
DeffsUnc = dat_mat(:,10);
Ds = dat_mat(:,11);
DsUnc = dat_mat(:,12);
nus = dat_mat(:,13);
nusUnc = dat_mat(:,14);

%hydro parametrization #2 (Gamma, Dmicro, Deff, nu)
Gammas2 = dat_mat(:,15);
Gammas2Unc = dat_mat(:,16);
Deffs2 = dat_mat(:,17);
Deffs2Unc = dat_mat(:,18);
Ds2 = dat_mat(:,19);
Ds2Unc = dat_mat(:,20);
nus2 = dat_mat(:,21);
nus2Unc = dat_mat(:,22);

%pheno parametrization (A0, A1, B0, B1)
A0s = dat_mat(:,23);
A0sUnc = dat_mat(:,24);
A1s = dat_mat(:,25);
A1sUnc = dat_mat(:,26);
B0s = dat_mat(:,27);
B0sUnc = dat_mat(:,28);
B1s = dat_mat(:,29);
B1sUnc = dat_mat(:,30);


%Plot results
FigHandle_Hydro = figure('name','Hydro Model');
NRows = 2;
NCols = 3;
subplot(NRows,NCols,1);
errorbar(Ts, Ds*1e3, DsUnc*1e3, DsUnc*1e3, TsUnc, TsUnc, 'r.');
% errorbar(Ts,Ds*1e3,DsUnc*1e3,'r.-');
xlabel('T (t)');
ylabel('D_{inc} (latt^2/ms)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('Incoherent/Microscopic D');

subplot(NRows,NCols,2);
errorbar(Ts, Gammas*1e3, GammasUnc*1e3, GammasUnc*1e3, TsUnc, TsUnc, 'r.');
% errorbar(Ts,Gammas*1e3,GammasUnc*1e3,'r.-');
xlabel('T (t)');
ylabel('\Gamma (KHz)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('Momentum Rel Rate');


subplot(NRows,NCols,3);
errorbar(Ts, nus*1e3, nusUnc*1e3, nusUnc*1e3, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('\nu (latt^2/ms)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('Viscosity');


subplot(NRows,NCols,4);
errorbar(Ts, Deffs*1e3, DeffsUnc*1e3, DeffsUnc*1e3, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('c (latt/ms)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('Deff');


subplot(NRows, NCols, 5);
% DFull = Ds.^2./Gammas + Deffs;
%I don't think this error propogation is really fair for D. These variables
%are correlated, and the error in DFull is actually quite a bit smaller.
%TODO: better errorbar for D
% csSqrtDivGammaUnc = Ds.^2./Gammas.*sqrt((DsUnc./Ds).^2 + (GammasUnc./Gammas).^2);
% DFullUnc = sqrt(csSqrtDivGammaUnc.^2 + DeffsUnc.^2);
errorbar(Ts, Deffs*1e3, DeffsUnc*1e3, DeffsUnc*1e3, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('D (latt^2/ms)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('D_{eff}');


subplot(NRows,NCols,6);
OneOverDUnc = 1./Deffs.^2.*DeffsUnc;
errorbar(Ts, 1./Deffs/1e3, OneOverDUnc/1e3, OneOverDUnc/1e3, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('1/D_{eff} ms/latt^2)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('1/D_{eff}');

suptitle(sprintf('Hydro Model Params Vs. T\n dt^2n + (\\Gamma + (\\nu + D)*k^2)dtn + (c^2 + \\Gamma*D + \\nu*D*k^2)*k^2n = 0'));

if save_results
    savefig(FigHandle_Hydro,'HydroModel_Summary.fig');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FigHandle_Hydro2 = figure('name','Hydro Model 2');
NRows = 2;
NCols = 3;
subplot(NRows,NCols,1);
errorbar(Ts, Ds2*1e3, Ds2Unc*1e3, Ds2Unc*1e3, TsUnc, TsUnc, 'r.');
% errorbar(Ts,Ds*1e3,DsUnc*1e3,'r.-');
xlabel('T (t)');
ylabel('D_{inc} (latt^2/ms)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('Incoherent/Microscopic D');

subplot(NRows,NCols,2);
errorbar(Ts, Gammas2*1e3, Gammas2Unc*1e3, Gammas2Unc*1e3, TsUnc, TsUnc, 'r.');
% errorbar(Ts,Gammas*1e3,GammasUnc*1e3,'r.-');
xlabel('T (t)');
ylabel('\Gamma (KHz)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('Momentum Rel Rate');


subplot(NRows,NCols,3);
errorbar(Ts, nus*1e3, nusUnc*1e3, nusUnc*1e3, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('\nu (latt^2/ms)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('Viscosity');


subplot(NRows, NCols, 4);
% errorbar(Ts, Ds*1e3, DsUnc*1e3, DsUnc*1e3, TsUnc, TsUnc, 'r.');
% xlabel('T (t)');
% ylabel('D_m (latt^2/ms)');
% grid on;
% ax = gca;
% ax.YLim(1) = 0;
% xlim([0, 10]);
% title('Dmicro');


subplot(NRows, NCols, 5);
errorbar(Ts, Deffs2*1e3, Deffs2Unc*1e3, Deffs2Unc*1e3, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('D (latt^2/ms)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('D_{eff}');


subplot(NRows, NCols, 6);
OneOverDUnc = 1./Deffs2.^2.*Deffs2Unc;
errorbar(Ts, 1./Deffs2/1e3, OneOverDUnc/1e3, OneOverDUnc/1e3, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('1/D_{eff} ms/latt^2)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('1/D_{eff}');

suptitle(sprintf('Hydro Model2 Params Vs. T\n dt^2n + (\\Gamma + (\\nu + D)*k^2)dtn + (\\Gamma*D_{eff} + \\nu*D*k^2)*k^2n = 0'));

if save_results
    savefig(FigHandle_Hydro, 'HydroModel2_Summary.fig');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

FigHandle_Pheno = figure('name', 'Pheno Model');
NRows = 2;
NCols = 3;
subplot(NRows, NCols, 1);
errorbar(Ts ,A0s*1e3, A0sUnc*1e3, A0sUnc*1e3, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('A0s (KHz)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('A0');

subplot(NRows, NCols, 2);
errorbar(Ts, A1s*1e3, A1sUnc*1e3, A1sUnc*1e3, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('A_1 (latt^2/ms)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0,10]);
title('A1s');


subplot(NRows,NCols,3);
errorbar(Ts, B0s*1e6, B0sUnc*1e6, B0sUnc*1e6, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('B0s (latt^2/ms^2)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0, 10]);
title('B0s');


subplot(NRows, NCols, 4);
errorbar(Ts, B1s*1e6, B1sUnc*1e6, B1sUnc*1e6, TsUnc, TsUnc, 'r.');
xlabel('T (t)');
ylabel('B1s (latt^4/ms^2)');
grid on;
ax = gca;
ax.YLim(1) = 0;
xlim([0, 10]);
title('B1s');

suptitle(sprintf('Pheno Model Params Vs. T\n dt^2n + (A_o + A_1*k^2)dtn + (B_o + B_1*k^2)*k^2n = 0'));

if save_results
    savefig(FigHandle_Pheno, 'PhenoModel_Summary.fig');
end